subroutine show_hierarchy(N_hier)
  use globals
  use particles
  implicit none

  integer*4 :: N_hier,ilow,ipart
  real*4, dimension(N_hier) :: Nlow
  real*4 :: small_val

  Nlow(:) = 0
  do ilow=1,N_hier
     small_val = 1.0e10
     if(ilow==1)then
        do ipart=1,NumIn0
           if(In0%Mass(ipart) < small_val) small_val = In0%Mass(ipart)
        enddo
        do ipart=1,NumIn1
           if(In1%Mass(ipart) < small_val) small_val = In1%Mass(ipart)
        enddo
        do ipart=1,NumIn2
           if(In2%Mass(ipart) < small_val) small_val = In2%Mass(ipart)
        enddo
        do ipart=1,NumIn3
           if(In3%Mass(ipart) < small_val) small_val = In3%Mass(ipart)
        enddo
        do ipart=1,NumIn4
           if(In4%Mass(ipart) < small_val) small_val = In4%Mass(ipart)
        enddo
        do ipart=1,NumIn5
           if(In5%Mass(ipart) < small_val) small_val = In5%Mass(ipart)
        enddo
        Nlow(ilow) = small_val
     else
        do ipart=1,NumIn0
           if(In0%Mass(ipart) < small_val .and. In0%Mass(ipart) > Nlow(ilow-1)) small_val = In0%Mass(ipart)
        enddo
        do ipart=1,NumIn1
           if(In1%Mass(ipart) < small_val .and. In1%Mass(ipart) > Nlow(ilow-1)) small_val = In1%Mass(ipart)
        enddo
        do ipart=1,NumIn2
           if(In2%Mass(ipart) < small_val .and. In2%Mass(ipart) > Nlow(ilow-1)) small_val = In2%Mass(ipart)
        enddo
        do ipart=1,NumIn3
           if(In3%Mass(ipart) < small_val .and. In3%Mass(ipart) > Nlow(ilow-1)) small_val = In3%Mass(ipart)
        enddo
        do ipart=1,NumIn4
           if(In4%Mass(ipart) < small_val .and. In4%Mass(ipart) > Nlow(ilow-1)) small_val = In4%Mass(ipart)
        enddo
        do ipart=1,NumIn5
           if(In5%Mass(ipart) < small_val .and. In5%Mass(ipart) > Nlow(ilow-1)) small_val = In5%Mass(ipart)
        enddo
        Nlow(ilow) = small_val
     endif
    enddo
  
  write(logunit,*)' -particle mass hierarchy:'
  do ilow=1,N_hier
     write(6,*)'  ',ilow,Nlow(ilow)
  enddo
  call flush(logunit)

end subroutine show_hierarchy


! ****************************************************************************


subroutine determine_new_arrays(cut_fac1,cut_fac2)
  use globals
  use particles
  implicit none



  integer*4 :: ipart
  real*4 :: cut_fac1,cut_fac2
  real*4 :: min_mass

  min_mass = minval(In1%Mass)
  NumOut0 = 0
  NumOut1 = 0
  NumOut2 = 0
  NumOut3 = 0
  NumOut4 = 0
  NumOut5 = 0
  do ipart=1,NumIn0
     if(In0%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In0%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In0%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  do ipart=1,NumIn1
     if(In1%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In1%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In1%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  do ipart=1,NumIn2
     if(In2%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In2%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In2%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  do ipart=1,NumIn3
     if(In3%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In3%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In3%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  do ipart=1,NumIn4
     if(In4%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In4%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In4%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  do ipart=1,NumIn5
     if(In5%Mass(ipart) < 1.01*min_mass) then
        NumOut1 = NumOut1 + 1
     else if(In5%Mass(ipart) < cut_fac1*min_mass) then
        NumOut2 = NumOut2 + 1
     else if(In5%Mass(ipart) < cut_fac2*min_mass) then
        NumOut3 = NumOut3 + 1
     else
        NumOut5 = NumOut5 + 1
     endif
  enddo
  
  write(logunit,*)' -new NumPart_Total:'
  write(logunit,*)'  ',NumOut0,NumOut1,NumOut2,NumOut3,NumOut4,NumOut5,&
       'total=',NumOut0+NumOut1+NumOut2+NumOut3+NumOut4+NumOut5
  call flush(logunit)

end subroutine determine_new_arrays


! ****************************************************************************


subroutine populate_new_arrays(cut_fac1,cut_fac2)
  use globals
  use particles
  implicit none

  integer*4 :: ipart
  real*4 :: cut_fac1,cut_fac2
  real*4 :: min_mass

  write(logunit,*)' -populating new arrays'
  call flush(logunit)

  min_mass = minval(In1%Mass)
  NumOut0 = 0
  NumOut1 = 0
  NumOut2 = 0
  NumOut3 = 0
  NumOut4 = 0
  NumOut5 = 0

  do ipart=1,NumIn0
     if(In0%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In0%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In0%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In0%Mass(ipart)

     else if(In0%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In0%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In0%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In0%Mass(ipart)

     else if(In0%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In0%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In0%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In0%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In0%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In0%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In0%Mass(ipart)

     endif
  enddo

  do ipart=1,NumIn1
     if(In1%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In1%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In1%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In1%Mass(ipart)

     else if(In1%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In1%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In1%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In1%Mass(ipart)

     else if(In1%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In1%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In1%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In1%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In1%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In1%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In1%Mass(ipart)

     endif
  enddo

  do ipart=1,NumIn2
     if(In2%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In2%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In2%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In2%Mass(ipart)

     else if(In2%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In2%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In2%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In2%Mass(ipart)

     else if(In2%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In2%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In2%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In2%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In2%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In2%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In2%Mass(ipart)

     endif
  enddo

  do ipart=1,NumIn3
     if(In3%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In3%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In3%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In3%Mass(ipart)

     else if(In3%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In3%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In3%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In3%Mass(ipart)

     else if(In3%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In3%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In3%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In3%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In3%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In3%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In3%Mass(ipart)

     endif
  enddo

  do ipart=1,NumIn4
     if(In4%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In4%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In4%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In4%Mass(ipart)

     else if(In4%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In4%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In4%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In4%Mass(ipart)

     else if(In4%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In4%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In4%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In4%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In4%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In4%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In4%Mass(ipart)

     endif
  enddo

  do ipart=1,NumIn5
     if(In5%Mass(ipart) < 1.01*min_mass) then

        NumOut1 = NumOut1 + 1
        Out1%Coordinates(:,NumOut1) = In5%Coordinates(:,ipart)
        Out1%Velocity(:,NumOut1)       = In5%Velocity(:,ipart)
        Out1%Mass(NumOut1)              = In5%Mass(ipart)

     else if(In5%Mass(ipart) < cut_fac1*min_mass) then

        NumOut2 = NumOut2 + 1
        Out2%Coordinates(:,NumOut2) = In5%Coordinates(:,ipart)
        Out2%Velocity(:,NumOut2)       = In5%Velocity(:,ipart)
        Out2%Mass(NumOut2)              = In5%Mass(ipart)

     else if(In5%Mass(ipart) < cut_fac2*min_mass) then

        NumOut3 = NumOut3 + 1
        Out3%Coordinates(:,NumOut3) = In5%Coordinates(:,ipart)
        Out3%Velocity(:,NumOut3)       = In5%Velocity(:,ipart)
        Out3%Mass(NumOut3)              = In5%Mass(ipart)

     else

        NumOut5 = NumOut5 + 1
        Out5%Coordinates(:,NumOut5) = In5%Coordinates(:,ipart)
        Out5%Velocity(:,NumOut5)       = In5%Velocity(:,ipart)
        Out5%Mass(NumOut5)              = In5%Mass(ipart)

     endif
  enddo

  ! Generate new ParticleIDs
  do ipart=1,NumOut1
     Out1%ParticleIDs(ipart) = ipart
  enddo
  do ipart=1,NumOut2
     Out2%ParticleIDs(ipart) = NumOut1 + ipart
  enddo
  do ipart=1,NumOut3
     Out3%ParticleIDs(ipart) = NumOut1 + NumOut2 + ipart
  enddo
  do ipart=1,NumOut5
     Out5%ParticleIDs(ipart) = NumOut1 + NumOut2 + NumOut3 + ipart
  enddo

  write(logunit,*)'  New mass ranges:'
  write(logunit,*)'  -PartType1:',minval(Out1%Mass),maxval(Out1%Mass)
  write(logunit,*)'  -PartType2:',minval(Out2%Mass),maxval(Out2%Mass)
  write(logunit,*)'  -PartType3:',minval(Out3%Mass),maxval(Out3%Mass)
  write(logunit,*)'  -PartType5:',minval(Out5%Mass),maxval(Out5%Mass)
  call flush(logunit)

  write(logunit,*)'  New ParticleIDs ranges:'
  write(logunit,*)'  -PartType1:',minval(Out1%ParticleIDs),maxval(Out1%ParticleIDs)
  write(logunit,*)'  -PartType2:',minval(Out2%ParticleIDs),maxval(Out2%ParticleIDs)
  write(logunit,*)'  -PartType3:',minval(Out3%ParticleIDs),maxval(Out3%ParticleIDs)
  write(logunit,*)'  -PartType5:',minval(Out5%ParticleIDs),maxval(Out5%ParticleIDs)
  call flush(logunit)
  

end subroutine populate_new_arrays





